Variational matrix product ansatz for dispersion relations 
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A variational ansatz for momentum eigenstates of translation invariant quantum spin chains is 
formulated. The matrix product state ansatz works directly in the thermodynamic limit and allows 
for an efficient implementation (cubic scaling in the bond dimension) of the variational principle. 
Unlike previous approaches, the ansatz includes topologically non-trivial states (kinks, domain walls) 
for systems with symmetry breaking. The method is benchmarked using the spin-J XXZ antiferro- 
magnet and the spin-1 Heisenberg antiferromagnet and we obtain surprisingly accurate results. 



The density matrix renormalization group (DMRG) 
has proven to be the most successful variational method 
for strongly correlated quantum lattice systems in one 
spatial dimensionP^The associated variational class is the 
class of matrix product states (MPS) pi which has been 
generalized to higher dimensional systems^ and can also 
be applied directly in the thermodynamic limit! 4 *' It is 
now understood that the success of these ansatze can 
be attributed to the fact that they occupy that corner 
of Hilbert space that is characterized by an area scaling 
law of entanglement entropy, a property also satisfied by 
groundstates of gapped short-ranged quantum systems! 5 -! 
The same argument applies equally — under some gen- 
eral constraints — to the lowest excited states of such 
systemsP 

Accurate information about the lowest-lying excited 
states is important to relate theoretical models to exper- 
imental measurements via spectral functions. Excited 
states appear as poles in these spectral functions, with 
corresponding residues given by the spectral weight. In 
dynamic DMRG methods,^ information about the spec- 
trum of excited states is gathered from an approximate 
computation of the spectral function. The latest state- 
of-the-art algorithms first generate time-dependent cor- 
relation functions through a dynamic real-time evolution, 
after which highly accurate information about the spec- 
trum can be extracted using spectral analysis P However, 
this approach is limited by the fact that only reasonably 
short time scales are computationally accessible, due to 
the linear growth of entanglement under real-time evo- 
lution. This results in a broadening of the exact poles 
in spectral functions. Extracting high quality informa- 
tion requires a combination of working with a large bond 
dimension D w 0(1O 3 ), linear prediction to extend the 
range of accessible time scales and complex statistical 
machinery to extract the precise position of the pole. 
However, since low-lying excited states also satisfy an 
area law for the scaling of entanglement entropy, it should 
be possible to construct a more direct and efficient ap- 
proximation. 



Nevertheless, MPS-inspired variational ansatze for ex- 
cited states are rare. Most interesting is the case of 
translation-invariant states, where the Hamiltonian is 
block-diagonal in the different momentum sectors. Rom- 
mer and Ostlund proposed a Bloch-like ansatz that al- 
lowed them to get an early estimate of the Haldane gap 
in the spin-1 Heisenberg antiferromagnet P by adding a 
virtual boundary operator Q acting in the D-dimensional 
auxiliary space to the MPS approximation of the ground 
state and making a momentum superposition thereof. 
This is closely related to the general strategy of Bijl, 
Feynman, and Cohen, who act with the Fourier trans- 
form of a local physical operator O on the ground state 
to create excitations^. This strategy is called the single- 
mode approximation in the context of spin systemsP^ 
and it has been applied to MPS irP^. Other ansatze in- 
clude the projected entangled momentum states^ and 
very recently the proposal by Pirvu et aZ.P^ in which 
a momentum superposition is taken of the ground state 
MPS in which at a single site the set of matrices A 3 are 
replaced by the set of matrices B s that is variationally op- 
timized. This ansatz contains and extends the Rommer 
and Ostlund ansatz (B s = QA S ) and the single-mode ap- 
proximation (B s — Y^t ( s \0\t) A* 1 ). All of these propos- 
als exploit translational invariance on a finite lattice with 
periodic boundary conditions, which unfortunately intro- 
duces finite-size effects and prevents them from reaching 
the computational efficiency [0(D 3 ) with D the bond di- 
mension of the MPS] of indirect methods on systems with 
open boundary conditions. 

This paper introduces a variational ansatz that allows 
us to describe excited states directly in the thermody- 
namic limit. This ansatz generalizes,^ but differs by not 
relying on periodic boundary conditions, which is of key 
importance for the formulation of a computationally effi- 
cient implementation. This also allows for the possibility 
of topologically non-trivial excited states, which are very 
important in systems with symmetry breaking,^ but for 
which few direct alternatives are available. 
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We now consider a one-dimensional lattice of ci-level 
quantum systems described by a local, translation- 
invariant Hamiltonian H = £„ eZ T^T"" , with f the 
translation operator that shifts the lattice over a sin- 
gle site, and h an operator that acts nontrivially only 
on sites zero and one (we restrict to nearest-neighbor 
Hamiltonians for the sake of simplicity). We approximate 
translation-invariant ground states of such Hamiltonians 
with infinite size uniform MPS (uMPS), given by 

mA))=v{ (n„ ez £t=i^")^is>, 

where |s) = |. . . sis 2 • • ■), A s (s = 1, 2, . . . , d), constitute 
a set of D x D complex matrices acting on a D-level aux- 
iliary system, and i>l and are two D-dimensional vec- 
tors living at ±oo. The MPS construction has a gauge 
invariance under the gauge transform A s i— > GA S G~ 1 , 
with G an invertible D x D matrix. While leaving the 
gauge unspecified, we do assume that the transfer matrix 
E^ — Yl s =i A s ® A s has precisely one eigenvalue 1 with 
corresponding left and right eigenvectors (l\ and |r) of 
length D 2 , to which we can associate D x D matrices I 
and r, respectively, by reshaping them. These two matri- 
ces are Hermitian and positive and are assumed to be full 
rank. We choose the normalization (Z|r) = Tr(ir) = 1. 
In addition, we assume that all other eigenvalues of E\ 
lie strictly within the unit circle, so the spectral radius 
of E\ — \r)(l\ is smaller than PS. Under these condi- 
tions, the boundary vectors and vr do not feature in 
normalized expectation values of local operators. 

Within the philosophy of Bijl, Feynman and Cohen, a 
typical elementary excitations of a local gapped Hamilto- 
nian can be interpreted as a momentum superposition of 
a localized disturbance of the ground state. We therefore 
define a variational ansatz for excitations as 

= E GiKnfn vi (• • • A s ~ 1 B s °A Sl • • • ) v R \s), 

where A s and A s represent the same (for a topologi- 
cal^ trivial excitation) or different (for a topologically 
non-trivial excitation in the case of symmetry breaking) 
set of matrices for which and |\1/(A)) are equally 

good {i.e. same energy) uMPS approximations of ground 
states of H {e.g. obtained using the imaginary time- 
dependent variational principle^). The state \<& K {B)) 
has momentum k 6 [— tt,tt). The set of D x D matrices 
B s (s = 1, 2, . . . , d) contains the only variational param- 
eters in our class; they are also denoted as a D 2 d vector 
B. All expectation values are quadratic in B, and we 
define 

($ K (B)\$ K ,(B')) = 2ttS{k - n')B^ K B', 
($ K {B)\H - H\$ K ,(B')) = 27t5(k - k')B^ k B' . 

For an infinite system size, the momentum eigenstates 
\<fr K {B)) cannot be normalized to 1 but rather satisfy 
a 5 normalization. This 5-function also appears in the 



expectation value of every translation invariant opera- 
tor. The energy expectation value ($ K (B)\H\§ K >{B')) 
has a contribution H ($ K (B)|$ K /(B')), where H = 
(V(A)\H\V(A)) = |Z| (V(A)\h\V(A)) is the diverging 
ground state energy, which was therefore subtracted in 
the definition of H«f . The spectrum of excitation en- 
ergies uj at momentum k can then be obtained by solv- 
ing the dD 2 -dimensional generalized eigenvalue problem 
(Hk*\ ISli*'). Since our variational space is a linear sub- 
space of the Hilbert space, this generalized eigenvalue 
system can be recognized as the Rayleigh-Ritz equation. 

However, it can easily be seen that \/X € C DxD , the 
choice B s = c iK A s X - XA S results in |$ K (B)) = 0. For 
any n ^ 0, or for n = and \^{A)) ^ |*(^)), our lin- 
ear parametrization B has D 2 linearly independent zero 
modes which can be eliminated by fixing a part of the 
variational parameters. For |<3>(A)) = j'i'(A)), there ex- 
ists a gauge transformation G such that A s = GA S G~ 1 
and one can see that X = G leads to B s — for k = 0. 
Hence, there are only D 2 — 1 linearly independent zero 
modes, but we can fix one additional variational param- 
eter in B by imposing (A)\$ K {B)) = 0. This orthog- 
onality constraint is automatically satisfied in all other 
cases. Hence, in all cases only {d — 1)D 2 variational pa- 
rameters remain. We prove elsewhere that this freedom 
in fixing some variational parameters is related to the 
gauge freedom in the original manifold of MPS, and that 
we can construct a linear representation B{x) in terms of 
a (d — 1)D X D matrix x containing the free variational 

parameters such that B (x^ B (y) — tr[x^yp^. In this 
study, we find the dD x (rf — l)D matrix Vl that contains 
an orthonormal basis for the null space of the D x dD ma- 
trix L with entries L a ^p s ^ = [{A s )H x / 2 ] a ^. Reshaping 
Vl such that [V£] a p = [Vl]( qs ), 7 (for all a = 1, . . . , D, 
s = 1, . . . , d, 7 = 1, . . . , {d — 1)-D), the representation is 
given by B s {x) = l~ 1 / 2 V£xr~ 1 / 2 , with I and r the left 

and right eigenvector of E\ and E^. We then obtain 
(seeSH) 

B{x)^H K B{y) = [(l\H*Mf\r) + (l\H™$\r) 

+ m B B { &{l-E A A )^Hii\r) 
+ {l\Hii{l-Ei)-'E^ x ]\r) 

+ e +iK {l\HH{l - Eir'Ei^l - e+™Ef)-iE B ^\r) 
+ e-^{l\HH{l ~ Eir'E^il - e-^Ei)^Ei (x) \f) 
+ e^{l\H^ {x) {l-e^Eir'Ef^\f) 
+ C -^{l\H^\l~e-^Ei)^E A [x) \f) 
+ ^{l\H^ x)A {l-e^E A A )-'Efy ) \f) 
+ e- 2 ^{l\H^ )A {l - e-^Ej)-iEj {x) \r)] , 
where E% = Y%=i As ® ^ and H cE 
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If.,t,u,v=i (s,t\h\u,v) (A"B«) ® (C D). When A ^A, 
which does not necessarily imply ^ \ty(A)), 

we can substitute A «— e lv A in order to obtain 
| = up to an infinite phase, which 

seems to indicate that the momentum label is completely 
arbitrary. This is an artifact of not having momentum 
in a system with open boundary conditions. It does not 
appear when A = A. This inconsistency can be solved 
by fixing ip such that the dominant eigenvalue (largest 
in magnitude) of E% is positive. When |*(j4)) = 
up to phase, the dominant eigenvalue of E$ is 1. By 
assumption, the dominant eigenvalue of E^ and are 
always 1. All inverses of (1 - E%), (1 - Ef), (1 - c iK E^) 

and (1 — c~ 1R E^) should be read as 'pseudo-inverses^ 
that act as zero in the eigenspace corresponding to the 
eigenvalue one of the transfer operator. In the case of 
symmetry breaking with ^ \^(A)), the spectral 

radius p(E^) < 1 and the expressions (1 — e 1 "^)" 1 
and (1 - e-'^E^)- 1 denote the full inverses. If A and 
A satisfy the properties that were outlined before, it is 
straightforward to prove that the corresponding uMPS 
approximate ground states with maximal symmetry 
breaking, i.e. they yield extremal values for the expec- 
tation value of the order parameter associated with the 
symmetry breaking. 

Excitation energies can thus be found from diagonaliz- 
ing the effective (d— l)D 2 x (d— 1)D 2 Hamiltonian defined 
with respect to entries of x and y in B(xyH K B(y) ) since 
the effective norm matrix is now the unit matrix. As 
shown in the results below, we are often interested in the 
lowest excitation energies. Using an iterative method for 
the different (pseudo)-inversions, the action of the effec- 
tive Hamiltonian can be implemented as an 0(D 3 ) oper- 
ation and can be combined with a sparse eigensolver. 

We now illustrate the power of our variational ansatz 
for excited states using some benchmark problems. The 
first Hamiltonian under consideration is the spin-j XXZ 
antiferromagnet in the symmetry-breaking phase A > 1. 
The ground states of maximal symmetry breaking are an- 
tifcrromagnetic and also break translational invariance. 
We therefore perform a spin-flip (<j x ) on every second 
site, in order to obtain 



-Hxxz = J ^ of of 



n u n+l 



j y a y 

J n a n+1 



The gap closes at the critical point A = 1, where we 
obtain the spin-J Heisenberg antiferromagnet. Fig. [T] dis- 
plays the full spectrum of excited states obtained with 
our ansatz at D — 33 (full diagonalization becomes com- 
putationally demanding for much larger values of D) for 
A = 4, resulting in (d — 1)D 2 = 1089 topologically triv- 
ial excitations (A s — A s ) and 1089 topologically non- 
trivial excitations (A s = J2t ( s \ aX \t) ^')- As pointed out 
irP^, the elementary particle excitations in the symmetry- 
broken phase are topologically non-trivial kinks, and all 
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FIG. 1. Spectrum of the lowest lying excitations of the spin-J 
XXZ antiferromagnet with anisotropy parameter A = 4 at 
D — 33. Black circles indicate topologically non-trivial exci- 
tations, gray squares indicate topologically trivial excitations. 
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FIG. 2. (a) Simulation results for as a function of 

the anisotropy A for various values of D ranging from 10 to 
400 (gray circles), as well as the exact result with the Bethe 
ansatz (black line), (b) Absolute error on the energy gap 
versus absolute error on the energy density for various values 
of A and D. 



topologically trivial excitations are compound states con- 
taining an even number of kinks. Not only do we recover 
the elementary kink, we also obtain a large set of points 
that fall within the two-particle (topologically trivial) 
and three-particle (topologically non-trivial) continuum. 
This happens because our ansatz contains a single per- 
turbation that can spread out over a region of 0(log d D) 
sites. The two-particle states are states with fixed total 
momentum K\ + K2 mod 2ir — k, but which consist of a 
superposition of relative momenta Ak = K2 — Ati so as to 
confine the two particles into the region allowed by the 
ansatz. 

We can assess the accuracy of our approach as a func- 
tion of D by comparing the energy gap A^ z (lowest 
excitation energy at k = or k = tt) with the exact 

value Axx Z . Because this gap belongs to a topologically 
non-trivial excitation that only comes in pairs on lattices 
with periodic boundary conditions, the value of the en- 
ergy gap calculated irP3 using the Bethe ansatz on a lat- 
tice with periodic boundary conditions is twice the exact 
value. As illustrated in Fig. [2] we can even obtain highly 
accurate values of the energy gap very close to the critical 
point, by going to larger values of the bond dimension D 
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using a sparse implementation. Note that errors on the 
elementary excitation are negative. This violation of the 
variational principle is caused by subtracting an estimate 
of the ground state energy H that is too large. Also note 
that the error on the gap scales as the square root of the 
error on the energy (density) for low values of D, but 
is proportional to this error for larger D, except at the 
critical point. 

Secondly, we study the spin-1 Heisenberg antiferro- 
magnet, which has a translation- invariant ground state in 
a symmetry-protected topological phase and is character- 
ized by the presence of the Haldane gappS There is a long 
history of numerical estim ates of the Haldane gap using 
a variety of methods^HH Yig. [3] shows the spectrum of 
topologically trivial excitations of the spin-1 Heisenberg 
antifcrromagnet, obtained using our ansatz for D = 30, 
where excitation energies are colored according to their 
degeneracy (which was obtained without using the sym- 
metry explicitly). As expected, the elementary excitation 
around momentum n = 7r is the S — 1 magnon excita- 
tion, and the corresponding energy at k = ir is the Hal- 
dane gap. This spectrum is in agreement with a previous 
proposal that was obtained using a dynamic DMRG sim- 
ulation on a finite chain of up to 400 sites with values for 
the bond dimension up to D = 2000 (second reference 
irP). The most accurate results for the Haldane gap are 
obtained using ground state DMRG on finite lattices in 
the second reference oP (A = 0.41047925(4) on a lattice 
of 400 sites with D = 500) and the last reference of" 
(A = 0.4104792485(4) on a lattice of 2048 sites with D up 
to 2700). With modest computational resources (solving 
a single eigenvalue problem itcratively) we obtain simi- 
lar results for various values of D up to D = 208, which 
were chosen so that all Schmidt values with a certain de- 
generacy (given by half-integral spin representations) are 
present, but without explicitly using the SU(2) symme- 
try. We can even improve the estimate for the Haldane 
gap by two more significant digits from a scaling anal- 
ysis in D: A = 0.410479248463^*}°! . Unlike in the 
second reference irP, where the rest of the dispersion re- 
lation was much less accurate than the gap, we can now 
expect roughly the same accuracy for all points of the 
dispersion relation where the elementary magnon exists. 
Around k = tt/4, the elementary magnon excitation is 
absorbed into the two-magnon continuum and becomes 



unstable against decay into two magnons. 

We have presented a variational algorithm, based on 
the matrix product state formalism, to determine topo- 
logically trivial and non-trivial excited states of one- 
dimensional quantum lattices, directly in the thermo- 
dynamic limit. We envisage that this set of excited 
states can also be used to accurately determine spec- 
tral functions, as requested irP2l. In addition, we expect 
that our proposal can be extended to the setting of two- 
dimensional lattice systems, by replacing a single tensor 
in the network of projected entangled-pair states^, or to 
the setting of quantum field theories, by building on the 
continuous MPS proposal. 
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FIG. 3. Spectrum of the lowest lying excitations of the spin-1 
Heisenberg antiferromagnet at D = 30, labeled by their spin 
S degeneracy (color online). 



ACKNOWLEDGMENTS 

J.H., H.V. and D.J.W. would like to thank R. Bertl- 
mann and F.V. for inviting them to the University of 
Vienna, where this research was initiated. D.J.W. is also 
grateful to J.H. and H.V. for their kind hospitality at the 
University of Ghent. Research supported by the Research 
Foundation Flanders (J.H.), the FWF doctoral program 
Complex Quantum Systems (W1210) (B.P.), the Science 
and Technology Facilities Council (D.J.W.), the DFG 
(FOR 635 and SFB 631) (J.I.C.), the EU Strep project 
QUEVADIS, the ERC grant QUERG and the FWF SFB 
grants FoQuS and ViCoM. 



1 S.R. White, Phys. Rev. Lett. 69 2863 (1992); 
U. Schollwock, Rev. Mod. Phys. 77 259 (2005). 

2 M. Fannes, B. Nachtergaele, R. F. Werner, Commun. 
Math. Phys. 144, 443 (1992); F. Verstraete, J. I. Cirac, 
V. Murg, Adv. Phys. 57, 143 (2008); J. I. Cirac, F. Ver- 
straete, J. Phys. A: Math. Theor. 42, 504004 (2009); 
U. Schollwock, Annals of Physics 326, 96 (2011). 

3 F. Verstraete and J.I. Cirac, arXiv:cond-mat/0407066 
(2004); V. Murg, F. Verstraete and J.I. Cirac, Phys. Rev. 



A 75, 033605 (2007). 

4 G. Vidal, Phys. Rev. Lett. 98, 070201 (2007); R. Orus and 
G. Vidal, Phys. Rev. B 78, 155117 (2008). 

5 F. Verstraete, J.I. Cirac, Phys. Rev. B 73, 094423 
(2006); T.J. Osborne, Phys. Rev. Lett. 97, 157202 (2006); 
M.B. Hastings, J. Stat. Mech. P08024 (2007); J. Eisert, 
M. Cramer, M.B. Plenio, Rev. Mod. Phys. 82, 277 (2010). 

6 L. Masanes, Phys. Rev. A 80, 052104 (2009) 



5 



7 K.A. Hallberg, Phys. Rev. B 52 R9827 (1995); S. Ra- 
masesha, S.K. Pati, H.R. Krishnamurthy, Z. Shuai and 
J.L. Bredas, Synthetic Met. 85, 1019 (1997); T.D. Kiihner 
and S.R. White, Phys. Rev. B 60, 335 (1999); E. Jeck- 
elmann, F. Gebhard and F.H.L. Essler, Phys. Rev. Lett. 
85, 3910 (2000); E. Jeckelmann, Phys. Rev. B 66, 045114 
(2002). 

8 R.G. Pereira, S.R. White and I. Affleck, Phys. Rev. Lett. 
100, 027206 (2008); S.R. White and I. Affleck, Phys. Rev. 
B77, 134437 (2008). 

9 S. Ostlund and S. Rommer, Phys. Rev. Lett. 75, 3537 
(1995); S. Rommer and S. Ostlund, Phys. Rev. B 55, 2164 
(1997). 

10 A. Bijl, J. de Boer and A. Michels, Physica 8, 655 (1941); 
R.P Feynman, Phys. Rev. 94, 262 (1954); R.P. Feynman 
and M. Cohen, Phys. Rev. 102, 1189 (1956). 

11 D.P. Arovas, A. Auerbach and F.D.M. Haldane, Phys. Rev. 
Lett. 60, 531 (1988). 

12 E. Bartcl, A. Schadschneider and J. Zittartz, EPJ B 31, 
209 (2003); S.G. Chung and L. Wang, Phys. Lett. A 373, 
2277 (2009). 

13 D. Porras, F. Verstraete, and J.I. Cirac, Phys. Rev. B 73, 
014410 (2006). 

14 B. Pirvu, J. Haegeman and F. Verstraete, Phys. Rev. B 
85, 035130 (2012). 

15 L.D. Faddeev and L.A. Takhtajan, Phys. Lett. A 85, 375 
(1981). 



States satisfying this condition are also known as pure 
finitely correlated states. 

J. Haegeman, J.I. Cirac, T.J. Osborne, I. Pizorn, H. Ver- 
schelde and F. Verstraete, Phys. Rev. Lett. 107, 070601 
(2011). 

In preparation. 

For k / 0, the operators (1 — e 1K E^) and (1 — e~ lK E^) are 
not really singular and could also be inverted fully. This 
will lead to the same result due to our representation of B. 
J. Des Cloizeaux and M. Gaudin, J. Math. Phys. 7, 1384 
(1966). 

F.D.M. Haldane, Phys. Lett. A 93, 464 (1983); Phys. Rev. 
Lett. 50, 1153 (1983). 

S.R. White and D.A. Huse, Phys. Rev. B 48, 3844 (1993); 

E. S. Sorensen and I. Affleck, Phys. Rev. Lett. 71, 1633 
(1993); M. Takahashi, Phys. Rev. B 48, 311 (1993); 
S. Todo and K. Kato, Phys. Rev. Lett. 87, 047203 (2001); 
H. Nakano and A. Terai, J. Phys. Soc. Jpn. 78, 014003 

(2009) ; H. Ueda and K. Kusakabe, Phys. Rev. B 84, 
054446 (2011). 

PE. Dargel, A. Wollert, A. Honecker, LP McCulloch, 
U. Schollwock, T. Pruschke, arXiv: 1203.2523. 

F. Verstraete and J.I. Cirac, Phys. Rev. Lett. 104, 190405 

(2010) ; T.J. Osborne, J. Eisert and F. Verstraete, Phys. 
Rev. Lett. 105, 260401 (2010). 



